
### Code to generate bootstrap.rda, which are bootstrapped MQ estimates

#setwd("~/Documents/working/dynIRT/mq")
set.seed(100)
library(dynIRT)

load("mqData2013.Rda")
load("mq_out.rda")
load("mq_in.rda")

votes <- mqdata[,1:45]

## remove unanimous cases
indic <- rep(TRUE, nrow(votes))
for (i in 1:nrow(votes)){
  meanval <- apply(votes[i,], 1, mean, na.rm=TRUE)
  if (meanval == 0 | meanval == 1){
    indic[i] <- FALSE
  }
}
votes <- votes[indic,]
mqdata <- mqdata[indic,]
vote_year <- mqdata$time
rc <- t(votes)


Ntrials <- 100
pseudo.rc <- vector("list", Ntrials)
vi.result <- vector("list", Ntrials)


for(trials in 1:100){

pseudo.rc[[trials]] <- matrix(NA, nrow=nrow(rc), ncol=ncol(rc))
idealpts <- lout$means$x
alpha.mq <- lout$means$alpha
beta.mq <- lout$means$beta

for(i in 1:nrow(pseudo.rc[[trials]])){
	for(j in 1:ncol(pseudo.rc[[trials]])){
		if(!is.na(rc[i,j])) pseudo.rc[[trials]][i,j] = 1*(pnorm(alpha.mq[j] + beta.mq[j]*idealpts[i,vote_year[j]]) > runif(1)) 
	}
}

## Re-remove unanimous votes that are newly unanimous from randomization
keepthese <- which(!(apply(pseudo.rc[[trials]],2,mean, na.rm=TRUE) %in% c(0,1)))

## Recode
pseudo.rc[[trials]][ pseudo.rc[[trials]] == 0] <- -1
pseudo.rc[[trials]][is.na(pseudo.rc[[trials]])] <- 0

pseudo.data.in <- list(rc = pseudo.rc[[trials]][,keepthese],
			startlegis = data.in$startlegis,
			endlegis = data.in$endlegis,
			bill.session = data.in$bill.session[keepthese,,drop=FALSE],
			T = data.in$T) 

pseudo.cur <- list(alpha = cur$alpha[keepthese,,drop=FALSE],
		beta = cur$beta[keepthese,,drop=FALSE],
		x = cur$x)

pseudo.priors <- priors


vi.result[[trials]] <- dynIRT(.data = pseudo.data.in,
                    .starts = pseudo.cur,
                    .priors = pseudo.priors,
                    .control = {list(
                        threads = 1,
                        verbose = TRUE,
                        thresh = 1e-8,
			maxit=300
                        )})

} #end for(trials in 1:1000)


save(vi.result,file="bootstrap.rda")





